module funcs
    implicit none
    integer,parameter :: N=500!节点数
    character(len=16),parameter :: fname='../show/data.bin'
    real,parameter :: RE=6376.0e3,&!地球半径
        B0=3.12e-5,&!地球磁偶极子强度
        PI=3.1415926,&!pi
        MU0=PI*4.0e-7,&!真空磁导率mu0
        RHO=1.6726219e-33!磁场密度
    real,save :: L,f,Q0,Theta(N),s(N),xi(N),h(N),step
contains
    function funcs_s(kai)
    !计算磁力线长度。起点从赤道开始
        implicit none
        real,parameter :: sq3=sqrt(3.0)
        real :: kai(),funcs_s(),t()
        t=sqrt(1.0+3.0*kai**2.0)
        funcs_s=L*RE/2.0*(kai*t+log(sq3*kai+t))
    end function
    function funcs_B(theta)
    !计算余纬为theta处的磁感应强度
        implicit none
        real :: theta,funcs_B
        real :: c,s
        c=cos(theta)
        s=sin(theta)
        funcs_B=B0/L**3.0*(sqrt(1.0+3.0*c**2.0)/s**6.0)
    end function
    function funcs_P(t)
    !计算参数P
        implicit none
        real :: t,funcs_P
        funcs_P=-cos(t)/sin(t)
    end function
    function funcs_Q(t)
    !计算参数Q
        implicit none
        real :: t,funcs_Q
        funcs_Q=Q0*sin(t)**14.0
    end function
    subroutine init(f0)
        implicit none
        integer :: i
        step=1
        f=f0
        L=6.6
        Q0=(2.0*PI*f*RE/B0)**2.0*RHO*MU0*L**8.0
        Theta(1)=asin(1.0/sqrt(L))
        Theta(N)=PI-Theta(1)
        forall(i=2,N-1) Theta(i)=Theta(1)+(i-1)*(Theta(N)-Theta(1))/(N-1)
        s=funcs_s(cos(Theta))
        h=L*RE*sin(Theta)**3.0
    end subroutine
    subroutine savedata
        implicit none
        open(2333,file=fname,form='unformatted',access='sequential')
        write(2333)N,s,xi
        close(2333)
    end subroutine
end module
